#!/usr/bin/env python

import sys
import optparse
from itertools import izip

from metrics import theor_probability

def read_distribution(fname, stat_level):
    ##
    ## {0.95: {(10,200): 6}}
    ##
    fobj = open(fname)
    header = fobj.readline()
    if header.split()[:3] != ['#', 'n1', 'n2']:
        raise Exception("Bad header: %s" % header.strip())
    quantiles = map(float, header.split()[3:])
    try:
        q_id = quantiles.index(stat_level)
    except ValueError:
        raise Exception("There is not distribution for given stat_level: %s" % stat_level)
    distrib_dict = {}
    for line in fobj:
        spl = map(int, line.split())
        distrib_dict[tuple(spl[:2])] = spl[2:][q_id]
    fobj.close()
    return distrib_dict


def main():
    parser = optparse.OptionParser(
        usage = '%prog -j JOINED_FILE -d DISTRIBUTION_FILE -g GENOMES_NUM',
    )
    parser.add_option(
        '-j',
        dest = 'joined_file',
    )
    parser.add_option(
        '-d',
        dest = 'distrib_file',
    )
    parser.add_option(
        '-g',
        dest = 'genomes_num',
        type = 'int',
    )
    parser.add_option(
        '-l',
        dest = 'stat_level',
        type = 'float',
    )
    opts, args = parser.parse_args()
    if not (opts.joined_file and opts.distrib_file and opts.genomes_num and opts.stat_level):
        print >> sys.stderr, parser.usage
        sys.exit(1)
#    distrib_dict = read_distribution(opts.distrib_file, opts.stat_level)
    genomes_num  = opts.genomes_num
    
    with open(opts.joined_file) as joined_file:
        for line in joined_file:
            gene_ids, ints = line.strip().split('\t')
            ints = map(int, ints.split(','))
            
            n1 = min(ints[0], ints[1])
            n2 = max(ints[0], ints[1])
            n12 = int(ints[2])
            
            p = theor_probability(genomes_num, n1, n2, n12) 
            if p < (1.0 - opts.stat_level):
#            if n1 > 1 and n2 > 1 and n12 > n1-2 and n12 > n2-2:
                unknown_gene_id, known_gene_id = gene_ids.split(',')
                sys.stdout.write("%s\t%s\t%.10f\n" % (unknown_gene_id, known_gene_id, p))


if __name__ == '__main__':
    import gc
    gc.disable()
    
    main()
